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The head-on coUision of two nonrotating axisymmetric equal mass black holes is treated numer- 
ically. We take as initial data the single parameter family of time-symmetric solutions discovered 
by Misner which consists of two Einstein-Rosen bridges that can be placed arbitrarily distant from 
one another. A number of problems associated with previous attempts to evolve these data sets 
have been overcome. In this article, we discuss our choices for coordinate systems, gauges, and the 
numerical algorithms that we have developed to evolve this system. 
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I. INTRODUCTION 

The coalescence of two black holes is considered to be one of the most promising sources of gravitational waves [|] . 
In a series of papers [^-^ we investigate numerically a special case of the black hole coalescence problem, namely the 
head-on collision of two equal mass black holes. Numerical computations are extremely difficult due to coordinate 
singularities, large gradients in the metric components, and numerical instabilities inherent in the two black hole 
spacetimes. Building on the work of Dewitt, Cadez, Smarr, and E ppl ey (|-|^ (henceforth abbreviated as DCSE) and 
the more recent work involving distorted, single black holes ^^|l2|, many of these numerical problems have been 
overcome in the present work. 

Our work (like that of DCSE) is based on studying evolutions of the Misner initial data set , a single parameter 
family of time symmetric solutions that allows the initial separation between the two black holes and the total ADM 
mass of the system to be specified. Section || discusses this data set briefly. In section II] we present the basic 



equations, formalism, coordinate system and gauge considerations. Section |rV| details some of the various numerical 
problems that we encountered along with the computational methods we developed to overcome them, drawing 
parallels between our methods to those of DCSE when appropriate. Section ^ reviews the general numerical methods 
we use to integrate the discrete Einstein equations. We also present convergence studies testing the reliability and 
robustness of our code. 

Since this paper is devoted primarily to numerical methods and tests of our code, we refrain from presenting results 
such as extracted waveforms, energy radiated, horizon oscillations, etc., except where we discuss convergence results 
for our code. Instead we refer the reader to the series of related papers for a detailed analysis of the physics of 

colliding black holes. 



II. THE MISNER INITIAL DATA 



There are a number of ways to construct initial data representing two black holes. One of the simplest to work 
with is the single parameter family of analytic data derived by Misner for the Einstein-Rosen model of two 
asymptotically flat sheets joined by two throats. Detailed studies of the Einstein-Rosen bridge construction and 
variations of it were discussed by Misner Lindquist Brill and Lindquist and others. Data sets of this 
type were first investigated numerically by Hahn and Lindquist and later by DCSE. 

The spatial 3-metric for the Misner data is written in the following conformally flat form using cylindrical coordinates 

dP = {dp^ + dz^ + p^d(t)^) . (1) 
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The conformal factor defined by 



*Af = 1 

and 
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solves the Hamiltonian constraint (|^) with the proper isometry imposed between the upper and lower sheets. This 
data set is both axisymmetric and time symmetric {Kij = 0) and represents two equal mass black holes with zero 
rotation. The two black hole centers are aligned along the axis of symmetry (z-axis) so the physical interaction is a 
head-on collision along the axis. The free parameter n is related to the physical parameters M (half the total ADM 
mass and approximately the mass of a single black hole when the holes are "infinitely" separated) 

oo 
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and L (the proper distance along the spacelike geodesic connecting the throats) 
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The effect of increasing is to set the two black holes further away from one another and decrease the total mass of 
the system. 

Comparing with the isotropic form of the Schwarzschild metric 
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we can see how the Misner data (|l|) is similar to (||). If one associates m\ ^ ^ 2/sinhn/i as the mass of a single 
black hole when the two holes are separated by large distances (/Lt —>■ oo), then the Misner data resembles (||) if we 
identify m = mi for regions close to one of the throats and m = 2mi in the far field. (See also, for example, Brill and 
Lindquist Q) 

On the initial time-symmetric surface, the throats are minimal area surfaces. As the initial data parameter fi 
is varied, the shape of the initial apparent horizon varies. If the holes are close enough together (small /i), a new 
minimal surface appears, surrounding both black holes on the initial slice. Cadez ||l^ has calculated the critical value 
fic — 1.362 when this occurs. We know that event horizons lie outside of, or are coincident with the outermost trapped 
surface ||2^. This implies that for values of fj, less than or equal to /ic the initial data is a single distorted black hole. 
As we discuss in [p|,p|pH| we have determined the fj, required for a single connected event horizon by integrating 
photons through the spacetimes. We find that for values of fj, greater than about 1.8 the black holes do not have a 
common event horizon on the initial timeslice. 



III. THEORETICAL FRAMEWORK 

We use the 34-1 (or ADM) formalism to write the general 4- metric as 

ds^ = -(a^ - f3'(3i)dt^ + 2l3Mdt -f jijdx'dx^, (7) 

where a is the lapse function foliating the four dimensional spacetime with three dimensional spatial hypersurfaces 
7ij, and /?' is the shift vector that specifies three dimensional coordinate transformations from time slice to time slice. 
Throughout this work we use Greek indices (ranging from to 3) to label four dimensional coordinates, and Latin 
indices (ranging from 1 to 3) to label spatial coordinates. We use geometric units in which the gravitational constant 
and the speed of light are set to unity. 

In the 3+1 formalism, the vacuum Einstein equations reduce to four constraint equations 

R - K,jK'^ = 0, (8) 



2 



Vj{K'^ -Y^K) = 0, 



(9) 



and twelve evolution equations 

9t7y = -2aXy + V,/3j + V^A, (10) 

+/3'=Vfci^y- + Kk.Vrp'' + K.kVjP^. (11) 

Here R is the Ricci scalar formed from the spatial 3-metric 7^ , Kij is the extrinsic curvature, K is the trace of Kij 
and Vj is the covariant derivative with respect to 7^^ . 



A. The Coordinate Systems 

Because the spacetimes we work with possess an axial Killing vector (which we set to be d/dx^), all variables are 
independent of x^. Denoting the azimuthal angle hy x^ = (p and using the standard {z,p,if)) cylindrical coordinates, 
we write the 3-metric for a general axisymmetric spacetime as 




7»j { b ) (12) 

.0 p^d, 

in the coordinate order a:' = (2, p, (f>). The variables a, b, c and d are functions of the coordinates z, p and t and are 
assumed to be asymptotically flat. The conformal factor ^ is a function of z and p only and does not evolve in time. 
It is determined on the initial time slice to satisfy the Hamiltonian constraint (|8|) and for the type of initial data that 
we consider, and all its derivatives are known analytically, given by equations (H) and (^. At the initial time slice, 
the 3-metric (|l^) takes the form of the time symmetric Misner data set (^ with a = h = d=l, c = Q and = m- 
The Einstein equations are simplified when a conformal factor is introduced into the extrinsic curvature in a manner 
similiar to the 3-metric (|l2|). We define the following form for the extrinsic curvature 

(ha K \ 

K,j=-^^k,,^-^^ \K hb ) . (13) 




The evolution equations ( |10| ) and ( |ll| ) can now be formulated as evolution equations for the metric components (a, 
6, c, d) and their corresponding curvature components (/la, hf,, he, hd). 

As we discussed in section the initial data we consider consists of two throats connecting two isometric sheets. 
By construction the throats are spheres, centered on the axis along p = 0, on which boundary conditions relating the 
metric across the two sheets may be imposed. Since the natural boundaries (the throats and a sphere surrounding the 
system far from the throats) do not lie along constant {z,p) coordinates, it is useful to introduce the "quasi-spherical" 
Cadez |5j coordinates {ri,£,) with -q being a logarithmic "radial" coordinate and ^ an "angular" coordinate. Cadez 
coordinates are related to cylindrical coordinates through the complex transformation 



X(C) =Tl + iS, 
_ 1 
~ 2 



i[ln(C + Co)+ln(C-Co)] (14) 
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where C, = z ~\- ip^ and Co ~ cothp is the value of C at the throat center. The constant 77 and C coordinate lines of 
( p^ ) lie along the field and equipotential lines of two equally charged metallic cylinders located at the centers of the 
two throats z = ± coth p. The coefficients C„ are determined by a least-squares method to set the throats (defined 
by -I- {zth ± coth/i)^ — 1/sinh^/i) to lie on an 77 = 770 = constant coordinate line. Both 770 and the different C„ 
are computed on a problem specific basis for different p using this least squares procedure. As the Cn are rapidly 
converging, the series ( p^ ) can be truncated to low order. In our simulations we typically keep only terms up to 
n - 15. 

The constant Cadez coordinate lines as viewed in the cylindrical coordinate system are shown in Fig. |l^. The grid 
in the Cadez coordinates is shown in the accompanying Fig. |l|b. The advantage afforded by this set of coordinates 
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is that they are spherical near the throats of the black holes and also far away in the wave zone, thus allowing us to 
deal with throat boundaries and asymptotic wave form extractions in a convenient way. Also, because the natural 
boundaries of the spacetime now lie on constant (77, ^) lines, our spacetimes are constructed and evolved in a manner 
similar to the single distorted black hole spacetimes in references . The disadvantage is that the transformation 

( p^ ) introduces a singular saddle point at the origin (z = p = 0) that is not pres ent in cylindrical coordinates. This 
creates certain numerical difficulties that we discuss in more detail in section 
The transformation satisfies the following Cauchy-Riemann conditions 

ri.p^-i.z-, n.zp = -i^zz ^ i,pp- (15) 

We note that these first and second derivatives of 77 and ^ can be computed "analytically" (to machine precision) from 
(p^). The Jacobian of the two coordinate systems after the Cauchy-Riemann conditions have been applied becomes 

J = {dT]/dpf + {dr]/dzf . (16) 



In this set of coordinates, we define the analog of the cylindrical based 3-metric (12) as 



(AC \ 

7*j- [C B , (17) 



sin^ ^ D , 



with the corresponding extrinsic curvature 



(Ha He \ 
K,, = He Hb , (18) 

V sin^ ^ HdJ 

in the order [r], ^, (f>]. The Misner initial data (0) can be written in the form of (|l7| ) using Cadez coordinates as 

ds^ = (drj^ + + sin^ ^^If—dcjjA , (19) 
V sin ^ / 

and identifying A = B = 1, C = 0, = JpV sin^ ^ and * = -i^M/J^/^. 

The success of our methodology depends critically on using both sets of coordinate systems (cylindrical and Cadez) 



to advantage as we discuss in section IV 



B. The Lapse and Shift 

Kinematic conditions for the lapse function a and shift vector complete the set of Einstein equations (^) through 
(prj). Even though the 3-metric is fixed on the initial slice by the Hamiltonian constraint, the lapse and shift can be 
chosen arbitrarily on the initial slice and thereafter. 

We impose the maximal slicing condition 

K^dtK = Q (20) 

throughout the evolution. Taking the trace of Eq. (|ll|) and inserting Eq. (^0|) results in the following elliptic equation 
for a 

VVia = aK,,K'\ (21) 

where we have used the Hamiltonian constraint to replace R with KijK^^. We choose to solve the lapse equation 
in the nonsingular cylindrical coordinate system for improved accuracy over Cadez coordinates. Solutions to the 
lapse equation written in cylindrical coordinates tend to be smoother and better behaved near the saddle point than 
solutions obtained by solving the equation written in Cadez coordinates. Since the evolution equations for the extrinsic 
curvature components contain second derivatives of the lapse, smoothness in this sensitive region is very important. 

For the initial lapse we have tried both a = 1, whereby observers are initially freely falling, and the solution of 
Cadez 11 
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+ ^ + - , (22) 

n=l ^ ' 

which is a generahzation of the standard static Schwarschild sUcing for a single black hole, but here the system is not 
static. In the first case {a = 1) the lapse is symmetric across the throat. In the second case it is antisymmetric and 
hence equal to zero on the throat. We find the antisymmetric lapse tends to work better as it "freezes" the evolution 
at the throat and slows it down in regions near and between the throats where calculations can be troublesome. For 
this reason, we work exclusively with antisymmetric boundary conditions for the lapse across the throat. 
The shift vector /3* is defined by imposing 

C^dtC = (23) 

to maintain a diagonal 3-mctric in the Cadez coordinates. This condition applied to Eq. ( p^ results in the differential 
equation 

Bi^+A^ = 2aHc, (24) 

for the nonvanishing shift vector components (3^ and /3^. We can rewrite Eq. (|2j) by introducing a shift potential fl 
that satisfies 



to get a single elliptic equation for 



We emphasize that it is the off-diagonal Cadez metric component C that is required to vanish and not the off-diagonal 
cylindrical component c. In general, the cylindrical metric is nondiagonal. This choice of gauge for the shift vector 
proved to be critical to the overall stability of the numerical evolution, particularly in suppressing the axis instability. 



IV. "PATCHED" COORDINATES 



We have investigated a number of different numerical schemes to solve the problem of colliding two black holes 
head-on. The basic idea that evolved from our progression of trials is to solve for the Cadez metric and extrinsic 
curvature components, defined by Eqs. ( p7| ) and (p^), on the Cadez grid and set C = dtC = 0. This approach has 
the advantage that the 3-metric is diagonal which helps to suppress the axis instability and simplifies the equations 
of evolution and the extraction of invariant gravitational waves in the far field. 

This approach has the further advantage that it is possible to define variables for the two black hole system that 
obey the same evolution equations with similar boundary conditions, as the single distorted black hole code developed 
in previous work ||l^,|ll|. In fact, the two black hole code in its final incarnation evolved from the code we developed 
for distorted axisymmetric single black hole spacetimes and much of the discussion in [ p^|JTl|] is directly applicable to 
the two black hole code. In this section we concentrate on those techniques developed and modified specifically for 
the two black hole problem and in particular the methods we use to overcome difficulties associated with the singular 
saddle point present in the Cadez coordinates. Further details relevant to both codes may be found in pO|,pT[|. 



A. The Grid 



Our spacetimes are equatorially plane symmetric, axisymmetric and isometric through a throat boundary. The 
computational grid is covered with Cadez coordinates and is bounded by the equator (z = 0), the axis (p = 0) and 
the isometry surface (r] = rjo). The outer boundary is set at 77 = rimax defined by ijmax — Vo = 5.8, which corresponds 
to an equivalent Schwarzschild radius ranging from ~ 125M to ~ 425M for ^ — 1.2 and /i — 3.25 respectively. The 
outer boundary is sufficiently far from the two throats that we can safely impose static conformal flatness boundary 
conditions on the metric and extrinsic curvature components in the asymptotic far field. 
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Other boundary conditions are specified in an analogous way to the single black hole spacetimes described in 
references Specifically, we require for all the Cadez grid based variables (except C, He and /3^) to be 

symmetric across the axis and equator. C, He and (3^ are antisymmetric across both the axis and equator. The 
throat isometry takes the form 

drjA ^ dr/B = d^rjD ^ dr,Hc = dr,P^ = (27) 

and 

a ^ C = Ha ^ Hb = Hd = f3'^ = (28) 

at ?7 = rjQ. 

These choices result from establishing a consistent convention to satisfy the Einstein equations subject to the 
symmetries and isometrics of the problem. The choice of lapse ( ^ ) for the initial data removes the freedom available 
to choose an isometry sign for the metric components and the boundary symmetric constraints on the shift help 
preserve the coordinate positions of the axis, equator and throat boundaries [[lOyil| ]. Boundary values for variables 
defined in the cylindrical coordinate basis are obtained from the tensor transformation of the corresponding Cadez 
grid variables. 

We make one final comment regarding the grid: the grid is comoving with the black holes throughout the entire 
evolution. The dynamical history is carried solely by the metric components. As the black holes approach each other 
the grid cells do not track the black holes and squeeze together in the impact area between the two throats. Instead 
the conformal metric function 7^^ — A, for example, goes to zero, signifying that the proper distance between grid 
lines is decreasing. 



B. The Saddle Point 



The most difficult problem associated with the Cadez coordinates is the coordinate singularity at the origin (z = 
p = 0), clearly evident in Fig. The difficulty with this point, which is inside the computational domain, is two 
fold: (?) The transformation Jacobian (16) between the Cadez coordinates {rj, ^) and cylindrical coordinates (z,p) 
goes to zero at the saddle point as (z^ + p^)Ko, where Ko is a constant defined in the neighborhood of the saddle 
point. Note that the Jacobian between cylindrical coordinates and Cartesian coordinates goes to zero just as ~ p at 
the saddle point, so that Cadez coordinates are singular with respect to Cartesian coordinates at z = p = 0. {ii) The 
77 axis turns abruptly with respect to the local three geometry at the saddle point. The 77 axis parallels the z axis for 
rj < rjg while becoming the p axis for rj > r]s, where rj^ is the 77 coordinate passing through the origin. 

The problem (i) associated with the vanishing of the Jacobian (volume element) is a familiar one. Similar problems 
exist in spherical and cylindrical coordinates at the origin. The problem is often treated by rescaling variables, 
factoring out the Jacobian in the numerical evolution. The second problem (ii) is unfamiliar and more troublesome. 
As the spatial three geometry is symmetric to reflection about the z = plane, and to rotation about the p axis, the 
time development of the geometries along the p axis and the z axis are intrinsically different. The holes are falling 
towards each other in the z direction. Although the initial data is smooth along the 77 axis, as soon as the holes 
start falling towards each other, the discontinuity in the metric functions, say jr^ri = A, will develop at 77 = 77s. For 
rj < rjs, jrjri decreases as the proper distance between grid points decreases. (In fact, there are two competing effects: 
the decrease in the proper distance between the holes, and the grid stretching effects whereby grid points along the 
77 axis with 77 < 77s fall towards the "north" hole.) On the other hand, 7^^ increases for 77 > 77s as grid points are 
falling towards the saddle point, again due to the grid stretching effect. The anisotropy in the p and z directions at 
the saddle point translates into the discontinuity along the 77 axis at rjs. The functions 7^^ on the two sides of rjg 
are really two different geometric objects. The spatial derivative driJnri is undefined analytically at 77 = 77s on the 
77 axis, and is very large in the finite differencing approximation. Likewise the finite differencing [jnriijl = rjs + Arj, 
^ « ^) — = '7s — At;, ^ « ^)]/(2A7;) for grid points off the i] axis but near the saddle point will also be large. 

To see this problem explicitly, we plot in Fig. ^ the extrinsic curvature function T^T^^ = Ha — — (9t7^^)/a evolved 
for just one time step. In this case, since the initial data and all spatial derivatives are known analytically, the 
extrinsic curvature can be computed analytically on the first time step (up to discretization in time), which can then 
be used to compute the metric function on the second time step without computing any numerical spatial derivatives. 
In this sense, the metric is known "analytically" , up to a finite difference in time. But the effect of the discontinuity 
at the origin (77 = 77s ~ 0.05, plotted as a circle in Fig. ||) is clear immediately. It is easy to see that the same 
difficulty exists for other quantities A, B , Hb - - etc. Without special precautions or treatments, the evolution quickly 
contaminates with numerical noise. It is possible that certain gauge conditions, such as quasi-isotropic [ p3[ or minimal 
distortion p^ , could minimize the effect of these discontinuities, but we have not tried them. The method we use in 
dealing with this problem is discussed in the following section. 
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C. Numerical Issues 



The most critical numerical issue is the treatment of the saddle point and the region near it. To help reduce this 
problem somewhat, we construct the Cadez coordinates in such a way as to avoid placing a grid point on the origin. 
Lines of constant ^ are staggered across the axis of symmetry and the equator so that both axis and equator lie on 
the half zones. Lines of constant rj are set relative to the throat position 77 = 770 at constant discrete spacing A77 
chosen to maintain a grid aspect ratio that is nearly unity. We do not enforce a constant 77 line to intersect the saddle 
point. Hence, the origin is straddled by both the angular and radial grid lines (see Fig. |l|a). Although this procedure 
eliminates some spurious effects arising from the singular point, it does not resolve the problem completely. Gradients 
in this vicinity are extremely large, for the reasons discussed above. 

The basic approach we use to evolve data near the saddle point is to take advantage of the fact that the cylindrical 
spacetime metric components ( |l^ ) are smooth everywhere, including the saddle point (although the volume element 
in cylindrical coordinates, or equivalently the Jacobian to Cartesian coordinates, is still zero at that point and along 
the z axis). We can therefore define a cylindrical coordinate "patch" to cover regions near the saddle point in the 
singular Cadez system. The patch is constructed beginning at the angular coordinate value ^ = 7r/2 and extended 
along the angular direction towards f = for a number of zones, depending on the angular resolution. The patch is 
also extended along the radial direction from the throat all the way out to the outer boundary to minimize distortions 
that might be suffered by radially propagating structures if they encounter patch boundaries "head-on" . In this patch 
region we evolve the cylindrical coordinate based metric (12) and extrinsic curvature components ( p^ ) on the Cadez 
grid and transform the solutions via the general tensor relations T[j = {dx^ / dx'^){dx^ / dx'^)Tu to reconstruct the 
Cadez components. 

We stress that this scheme is not a coordinate patch in the formal sense. Grid lines are no where laid along the 
(z,/)) coordinates, and derivatives are not taken in {z^p) coordinates. Rather, we are simply evolving two sets of 
components, Cadez and cylindrical, independently of each other (except for the coupling at the patch boundaries) on 
a single Cadez grid. The nonsingular cylindrical components are used to correct the singular Cadez components in 
the patched region. On the other hand, the Cadez components provide corrections to their cylindrical counterparts 
everywhere else, helping to suppress the axis instability that is inherently present in the cylindrical coordinate system 
possessing a nondiagonal metric. To help integrate the patched region into the rest of the spacetime for smooth 
evolutions, we construct a layer of buffer zones surrounding the patched region. Within this boundary of zones, both 
sets of components are evolved and a linear weighting scheme is used to blend all evolved variables to the values at 
the edges of this buffer domain. 

In a typical calculation, the lapse collapses approximately "spherically" along constant 77 lines, thereby freezing 
the fields interior to this domain. (This is the generic behavior of the maximal slicing condition or any singularity 
avoiding lapse function.) Hence, it is necessary to evolve with the patch in place only until the lapse drops below 
a critical value (typically 0.025) at the origin to prevent any evolution from ocurring there. Once the lapse falls 
below this value, the simulation is continued by evolving only the Cadez variable components over the entire Cadez 
grid, including the saddle area. The removal of the patch at relatively late times is necessary to maintain stable 
and accurate evolutions. Evolving the natural Cadez metric components on the Cadez grid does not suffer from the 
numerical instabilities inherited from applications of chain rule derivatives in regions of extreme gradients. 

We end this section with a discussion of an alternative strategy that we tried for the numerical evolution. Although 
it was unsuccessful, it is instructive to point out why it failed. Given that the cylindrical metric variables are well 
behaved near the saddle point, one might consider evolving the entire system in those variables. It is well known that 
the axis instability can be suppressed by choosing the shift vector (/?'', 0) such that c = 0. This strategy is effective 
at minimizing problems related to both the axis instability and the saddle point, but it introduces a new problem 
related to grid stretching. As noted above, very large gradients develop in the radial metric function surrounding 
the hole. This radial metric function is composed of the three cylindrical metric functions ^zz, Ipzi and that are 
actually being evolved. Using a shift that forces c = in this coordinate system causes extreme angular gradients to 
develop near the transition between where is primarily radial (along the 2;-axis) and where 7pp is primarily radial 
(along the equator). Therefore, along a line of about 45 degrees between the axis and equator, instabilities develop 
as the grid stretching becomes severe and the metric functions 7^2 and 7pp develop stepfunction-like discontinuities. 
A nonvanishing ^pz component serves to absorb some of this shear but introduces instabilities on the axis. For this 
reason, we adopted the hybrid scheme described above, whereby the Cadez metric variables are evolved over much 
of the grid, and the cylindrical metric variables are evolved over a smaller region covering the saddle. In this way 
we were able to benefit from the advantages afforded by each coordinate system, while minimizing the problems that 
each presents. 
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V. NUMERICAL METHODS AND CODE TESTS 



A. Solving the Discrete Einstein Equations 



The numerical integration of the evolution equations ( |lOD and (11) is performed in an unconstrained manner 
because of practical computational time limitations. We do not enforce either the Hamiltonian (^) or the momentum 
constraints during the course of evolution except at the initial time slice. 

Eqs. (^0|) and (^l|) are solved using the standard time explicit second order accurate leap frog method whereby the 
extrinsic curvature components are staggered by a half step in time relative to the metric components. Schematically 
we have 

and 

= K^l + f + At 

i t, L \ L LI 

+ f/3rVi^r + K'VP'^ - VVa"+'/') Ai, (30) 



where subscripts i (superscripts n) refer to discrete spatial (temporal) positions. To maintain second order accuracy, 
variables are extrapolated to the proper time levels as needed using the second order formula 

^..+i/2^3^„_l^„-i^ (31) 

This method propagates gravitational waves with less dissipation and dispersion than other methods we have tried pl| . 

Explicit methods require stringent restrictions on the size of timesteps to maintain stability. A condition we found 
to provide a good balance between computational speed and accuracy is At = 4MA77, where M is half the ADM 
mass defined by (Q) and A77 is the spacing interval in the radial direction. 

As an added measure of stability, we introduce numerical diffusive terms to the discrete evolution equations. This 
effectively adds to the right-hand-side a term of the form A:V^7y and kW^Kij to the evolution equations (|l^ and ( |Tl| ) 
respectively. The coefficient k is chosen as small as possible to minimize its effect on the accuracy of solutions while 
enhancing the stability of the time integration. We construct k dimensionally in the manner k ^ (Ax)'^ / {2At) to 
scale appropriately with the grid parameters. The proportionality constant is typically chosen to be of order < 0.05. 
We have verified that the addition of these diffusive terms has little effect on the extracted radiation waveforms for 
the dominant £ = 2 modes (but the higher frequency = 4 is more sensitive to this procedure, as discussed below) . 

Spatial derivatives appearing as source terms in the discrete time evolution equations are differenced using both 
standard second and fourth order center differences to approximate drj and on the uniform Cadez grid. We find that 
fourth order center differences provide more accurate solutions evident in calculations of apparent horizon masses and 
gravitational wave form extractions. However, because black hole simulations can develop large gradients, fourth order 
differences are generally less stable than using second order centered differences. (See also Refs. [ pl]j25| for discussions 
of the effect of second and fourth order spatial derivatives.) The results presented in this paper and in the series of 
companion papers [&-Hwere obtained by center differencing to fourth order all spatial Cadez derivatives appearing 
in the evolution Eqs. ( |l0| ) and ([Tl|). Within the coordinate patch, it is necessary to compute spatial derivatives with 
respect to z and p of the cylindrical coordinate based variables. Finite differences of such derivatives {dzQ, dpa,.. .etc.) 
are computed on the Cadez grid using the chain rule formulas {dp = ij^pdrj + Jy.z^^, ...etc). 

The elliptic equations ( pl| ) and ( p6| ) for the lapse and shift potential respectively are discretized using second order 
center differences. The resulting coupled algebraic equations give rise to large sparsely banded matrices which we 
solve using an iterative two dimensional multigrid algorithm developed by Steve Schaffer p9| . 



B. Convergence Studies 

There are limitations to the grid resolutions that can be achieved. First, if the angular spacing A^ is too small, the 
zones poised next to the axis can trigger the axis instability during the evolution. Second, the phenomenon of "grid 
sucking" in which the black hole absorbs coordinate lines throughout the evolution is enhanced with resolution. As 
the resolution is increased, the peaks corresponding to the increased proper distance between nodes on the grid near 
the horizon become more pronounced and gradients grow ever sharper as one obtains more accurate solutions. These 
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sharp features are difficult to resolve numerically and are ultimately responsible for developing errors at late times, 
causing the code to crash. With our fixed mesh we cannot afford the computer time to add finer zones arbitrarily to 
accurately resolve the peaks with sufficient radial resolution to suppress numerical instabilities for the duration of a 
simulation. The convergence studies presented in this article include 100 (27), 200 (35), and 300 (55) radial (angular) 
zones on a uniformly spaced grid. 

For each parameter run (different values of ii) we extract both the £ = 2 and £ ~ A waveforms at radii of 30, 40, 50, 60, 
and TOM. (Coordinate positions corresponding to physical distances in units of M are approximated from the initial 
data in the asymptotically spherical far field as r ~ — ^E'^/A/.) We use results at each of these radii to check 

the propagation of waves and the consistency of our energy calculations. Table | shows the energy radiated across the 
five detectors at different grid resolutions for the case /x = 2.2. The convergence rate is at least quadratic in the grid 
spacing for all detectors. More specifically, we find the convergence between the 200 and 300 radial zone evolutions 
to be in the range 0.2-0.7% for each individual detector, with the larger deviations occurring for the innermost and 
outermost detectors. The median deviation is ~ 0.4%. We also find a general agreement among the different detectors 
at the same resolution. For example, the rms deviation (taken over all detectors) from the average radiated energy 
is ^ 5% for the 200 radial zone evolution. A maximum deviation of ~ 15% occurs between the two detectors furthest 
from one another. Reasons for the larger discrepancies among the inner- and outer- most detectors are the following: 
(i) the detectors closer in have a greater difficulty in separating gravitational wave signals from other parts of the 
gravitational field, (ii) the detectors further out arc more susceptible to numerical resolution problems and artificially 
induced diffusion, and (iii) the finite numerical run times allow for more of the gravitational wave train to pass through 
the inner detectors. 

In Fig. ^ we show the £ = 2 waveform extracted at r = 40M for the three different spatial resolutions of 100 (27), 
200 (35), and 300 (55) radial (angular) zones for the case /i = 2.2. Up until t ^ 150Af the results are quite similar in 
all three cases. After that time the low resolution run develops some difficulty due to the poorly resolved peak in the 
metric functions and becomes unstable. We note that with the higher resolution simulations one obtains a slightly 
better fit to the known quasinormal waveform, especially at late times. The higher resolved waveforms suffer less 
dispersion and damping attributable to numerical effects. Prior to i ~ I25M, the waveforms agree to within ^3% for 
all resolutions. Furthermore, the diffusion and patch parameters are different for each run, so the £ = 2 waveforms 
are quite stable with respect changes in computational parameters. In fact for a fixed resolution, varying the patch 
parameters (such as the size of the patch, the time at which it is lifted and the numerical diffusion) the waveforms 
vary by no more than ^10%, and typically by less than a few per cent. 

Next we show results from the more difficult £ — A extraction in Fig. ^ for the same ^ = 2.2 case. The £ = A 
mode is clearly more sensitive to the computational parameters than the £ = 2 extraction in both the signal preceding 
the strong quasinormal ringing (beginning at i ~ 75Af) and in the amplitude of the ringing signal. Also for a 
fixed resolution, the amplitude of the £ — A waveforms varies by about a factor of two with large changes in the 
computational parameters, with the largest effect coming from the patch parameters. As the energy output varies 
quadratically with the wave amplitude, the energy carried by the £ — A modes is uncertain to about an order of 
magnitude. In spite of these effects, a fit of the two lowest £ = A quasinormal modes to the high resolution run, shown 
in Refs. ||,^, is quite good. 

The reasons for the £ = A extraction to be less accurate are clear. In the first place, the amplitude of the £ = A 
component is much smaller than that of £ = 2, hence harder to extract from the background noise level. Moreover 
the more complicated angular distribution of the £ = A component needs more angular zones to be fully resolved than 
have been used in these runs. Even though we are unable to determine with great certainty the absolute amplitude 
of the £ = A signal with our present code, it is clearly seen in the data and does match the expected quasinormal 
frequency. On the other hand, we are confident in the accuracy of the larger and more robust £ = 2 signals, which 
are not sensitive to computational details as shown in Fig. |[ These figures show both the strengths and limitations 
of the present code. 

VI. SUMMARY 

In this paper we have presented the methods we developed to evolve the head-on collision of two equal mass black 
holes initially at rest. This problem is a difficult one that has required a number of numerical strategies designed 
to handle large gradients, to suppress instabilities, and to deal with singular points in the coordinate system. The 
result of our work is a code that can accurately evolve the black hole collision problem for a range of initial data 
sets. By considering simulations run with different numerical parameters and at different resolutions and we have 
demonstrated that the £ = 2 waveforms and hence the total radiated energy calculations are accurate. The amplitude 
and early time behavior of the £ = A waveforms are more sensitive to numerics, although the essential features (the 
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wavelength and damping time) of the quasinormal ringing are clearly present and can be accurately resolved. In 
a series of companion papers 0^ we present detailed analyses of the physical results obtained using the methods 
outlined in this article. 

Our work represents a step towards solving the more general and dynamic problem of coalescing binary black holes. 
It is our long term goal to develop a multi-purpose three dimensional code capable of simulating dynamic fully general 
relativistic spacetimes containing multiple black holes of arbitrary mass, rotation and impact parameters. 
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FIG. 1. (a) The Cadez grid is constructed for the case fi = 2.2 and displayed in a single quadrant with cylindrical coordinates. 
The throats are centered on the axis at z = ± coth jj,. Lines of constant 77 concentrically surround the throat locally, and become 
spherical far from the holes, (b) The computational grid is shown in Cadez coordinates. 
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FIG. 2. The extrinsic curvature function Jf^^ = Ha is plotted along a ^ = constant line (~ 7r/2) after a single time step in 
evolution. The geometric discontinuity in this function is evident. 

FIG. 3. We show the £ = 2 waveform at various resolutions of 100 (27), 200 (35), and 300 (55) radial (angular) zones 
for the case /i = 2.2. Except for the low resolution at very late times, the waveforms agree to within less than 3% across all 
simulations. 

FIG. 4. We show the i = 4 waveform at various resolutions of 100 (27), 200 (35), and 300 (55) radial (angular) zones for the 

case /i = 2.2. The amplitude of the ringing modes for these simulations varies by about a factor of 2. depending upon the grid 
resolution and other computational parameters such as the size and duration of the numerical patch and the added numerical 
diffusion. In spite of the uncertainty in the amplitude of the signals, the i = 4 quasinormal mode is clearly present in all cases. 



detector 


100 radial zones 


200 radial zones 


300 radial zones 


30M 


7.032 xlO"" 


6.098 xlQ-" 


6.068 xlO""* 


40M 


6.052 xlO"" 


5.785 xlQ--* 


5.773 xlO"" 


50M 


5.710 xlO"" 


5.619 xlQ--* 


5.606 xlO^"* 


60M 


5.346 xlO"" 


5.476 xlQ-* 


5.461 xlO"" 


70M 


5.069 xlO""* 


5.274 xlO"** 


5.313 xlO"* 



TABLE I. Convergence study of the total radiated energy for the case n = 2.2. The energies are computed at the five 
detector locations for three different spatial resolutions. The convergence rate is at least quadratic in the grid spacing for all 
detectors and deviations between the 200 and 300 radial zone simulations are on the order of a few percent. 
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